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We investigate the bandwidth-controlled Mott transition in the Hubbard model on the 
checkerboard lattice at half filling using the path-integral renormalization group (PIRG) 
method. It is demonstrated that the system undergoes a first-order phase transition to the 
plaquette-singlet insulating phase at a finite Hubbard interaction. This conclusion is drawn via 
a detailed analysis of the spin and charge correlations around the phase transition point by 
means of the PIRG method aided with a new iteration scheme introduced in this paper. 

KEYWORDS: Mott transition, Hubbard model, checkerboard lattice, geometrical frustration, path integral 
renormalization group method (PIRG) 
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1. Introduction 

Strongly correlated electron systems with geometrical 
frustration have attracted much interest recently. A well- 
known example is the frustrated pyrochlore lattice, which 
is given by a three-dimensional corner-sharing network 
of tetrahedra. This family includes the transition-metal 
oxides L1V2O4 1 and TI2RU2O7, 2 ' 3 where heavy fermion 
behavior and the Mott transition without magnetic or- 
dering were observed at low temperatures. In these com- 
pounds, electron correlations on the frustrated lattice 
may be a source of intriguing properties at low temper- 
atures. These experimental findings have stimulated the 
intensive studies of the Hubbard model on the pyrochlore 
lattice and its two-dimensional (2D) analog called the 
checkerboard lattice (Fig. 1). 4 ~ 16 We will focus on the 
checkerboard lattice model in this paper. 




Fig. 1. (Color online) Checkerboard lattice. Boxes indicate the 
clusters used in this paper and the number labels the size of 
them. 



The pyrochlore- and checkerboard-lattice models with 
nearest-neighbor hopping have some peculiar properties 
at half filling; the Fermi level just touches a flat band 
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(Fig. 2), which invalidates standard perturbation calcu- 
lations 4-6 and thus makes theoretical investigations ex- 
tremely difficult. This in turn poses a challenging the- 
oretical problem in the Mott physics in the frustrated 
systems. In the early work on the checkerboard lat- 
tice model, it was claimed that there may be a metal- 
insulator (MI) transition at infinitesimally small interac- 
tion, which is followed by the second insulator-insulator 
transition. 5 ' 6 Later on, however, it was demonstrated 
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Fig. 2. (a) The band structure along the symmetry lines in the 
Brillouin zone (B.z.) and (b) the density of states in the nonin- 
teracting case on the checkerboard lattice (eq. (2)). 



that the paramagnetic metallic state is stabilized at least 
in the weak coupling region, 17 ' 18 although it remains still 
open how the paramagnetic metallic state competes with 
other states, and how it is connected to the plaquette 
valence-bond crystal state realized in the strong coupling 
limit. 19 

In this paper, we investigate electron correlations in 
the Hubbard model on the checkerboard lattice at half 
filling at absolute zero by means of the path-integral 
renormalization group (PIRG) method. We show that 
the system undergoes a single first-order phase transition 
to the plaquette-singlet Mott phase at a finite Hubbard 
interaction. To carry out the calculation, we propose a 
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new iteration scheme for the PIRG algorithm, which sub- 
stantially improves the numerical accuracy, and hence al- 
lows us to figure out the nature of the phase transition. 

The paper is organized as follows. In §2, we introduce 
the model Hamiltonian and briefly explain the PIRG 
method. The newly introduced iteration scheme is also 
mentioned in the section. We discuss the phase transi- 
tion in the checkerboard Hubbard model in §3. A brief 
summary is given in §4. 

2. Model and Method 

We start with the Hubbard model on the checkerboard 
lattice (see Fig. 1), 

U=-t ^ima^jm'a + U Y ^ m \Tl lmh (1) 

where Cj TOCT (ct ) is an annihilation (creation) operator 
of electron in the z-th unit cell with spin a and sub- 
lattice index m (=1,2), and fii ma = c\ m<y Ci ma . U is the 
Coulomb interaction and t is the transfer integral with 
the same couplings on vertical, horizontal and diagonal 
bonds of the checkerboard lattice. A unitary transfor- 
mation of the kinetic term of the Hamiltonian Ji k gives 
the diagonalized form H k = J2k,a,a £ a(k)a\ aa a ka <n with 
two eigenvalues for each fc, 

e (k)-f 2t for a = 1, 

a( ' \ -2t(l + cos k x + cos k y ) for a = 2, 1 ' 

where a represents the band index. We assume t > 0, 
hereafter. The dispersion relations and the correspond- 
ing density of states are depicted in Fig. 2. The topmost 
energy band is flat over the whole Brillouin zone (B.z.), 
while the lower one is dispersive with the characteris- 
tic band structure for the square lattice with nearest- 
neighbor hopping. It is known that the band struc- 
ture originates from highly frustrated geometry of the 
checkerboard lattice. 20 

As discussed in Refs., 4 ' 5 the perturbation expansion 
in U encounters divergence at third and higher orders, 
which is due to the presence of the flat band in the sys- 
tem. Therefore, powerful numerical techniques are nec- 
essary to discuss the phase transitions in the model. In 
the frustrated system, the quantum Monte Carlo simu- 
lations usually suffer from the minus sign problems, and 
the exact diagonalization may not be efficient to discuss 
the ground state properties in the thermodynamic limit. 
Although one can deal with larger systems by means of 
the variational Monte Carlo simulations, the obtained re- 
sults strongly depend on a trial wave function employed. 
In this way, it may not be easy to discuss the nature 
of the phase transitions in the system. To overcome the 
difficulty, we here make use of the PIRG method devel- 
oped by Imada et al. 21,22 which allows us to improve 
the ground state systematically. It has an advantage in 
treating frustrated electron systems with large clusters, 
in contrast to the other numerical methods. 

The PIRG treatment is in principle independent of an 
initial state and an iterative method employed. However, 
careful choices of them are important to converge the 
PIRG calculations within available computational time. 
To obtain an appropriate initial state, we make use of 
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the unrestricted Hartree-Fock (UHF) approximation. 23 
In the approximation, the site- and spin-dependent mean 
fields (himcr) are introduced and the interaction term Jiu 
in the Hamiltonian is then replaced by 

)(«iml)], (3) 

where the parameters {hi ma ) are determined by the self- 
consistent equations {n ima ) = ((j)o\n ima \4> ) . Here |</>o) 
is the ground state of the Hamiltonian 'Hj. + 7t!^ hf . We 
use the resulting wave function as an initial one for the 
PIRG simulations, and then take into account quantum 
fluctuations. Note that the interaction C7j„ t is not nec- 
essarily equal to the original U when the initial wave 
function is determined. By performing the PIRG itera- 
tion, the approximate ground state is described by the 
Slater basis states as \ip) = Yla=i c a\4>a), where c a is 
an amplitude of \<fi a ), and L is the dimension of the 
truncated Hilbcrt space. Furthermore, by using an en- 
ergy variance extrapolation scheme, 21,22 ' 24 we can de- 
duce physical quantities such as the ground state energy 
and the double occupancy. Here the energy variance is 
defined by [(H 2 ) - (H) 2 ]/(H) 2 . 

We now introduce a new iteration scheme which can 
substantially improve the convergence of PIRG algo- 
rithm for the frustrated systems. Note that the true 
ground state could be obtained in principle by acting 
the time evolution operator on the initial state as 

|^> = e-^o>, (4) 

with imaginary time (5 — > oo. To approach the ground 
state in the framework of the PIRG method, we divide 
the operator into those in a small imaginary time slice 
At. In terms of the Hubbard-Stratonovich transforma- 
tions, the operator is explicitly given as 

e -ArW = ^ e -ArW t /2 y {{s}) £ -Arn k /2 

W 

+ 0((Ar) 3 , (5) 
N 1 

V(M) = Jie a (" ) **Te«(-»*)*u ) (6) 
»=x 

where £ {s} represents £ Sl =±i E a ±1 ■ ■ ■ £ flJf =± n 
a(s) = 2as - AtU/2 and a = tantT 1 i/tanh (Ar[//4). 

In the ordinary PIRG method, the Hilbert sub- 
space with L basis states is generated by the time 
evolution operators for the local interacting term 
e - ATUn ir n n . Since the operator is replaced by the po- 
tential terms with one auxiliary field as e - ATUn ii n ii = 
1 / 2 E <( =±i ea( * , « a( "" l6il i i1: generates two basis 
states. Therefore, by acting the operators with differ- 
ent i on a certain basis state a couple dozen times, 
we can easily produce L states. However, it is difficult 
to numerically generate independent basis states in this 
procedure, because the effective dimension is often de- 
creased. This drawback may result from the fact that 
two states generated by the local Hubbard-Stratonovich 
field are very similar, making it sometimes difficult to 
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distinguish themselves within given numerical accuracy. 
This may cause the reduction in the effective dimension 
of the Hilbert subspace. To improve the situation, we 
here select L basis states randomly within the 2 N states 
generated by the operator e~ Ar H with the full Hamil- 
tonian. In this case, unfavorable history accumulated via 
successive actions of the projection operators does not 
appear and thereby we can effectively avoid the decrease 
in the dimension of Hilbert subspace. In particular, when 
At' 0) is small, the above procedure improves the 
PIRG results substantially, and the introduction of ran- 
dom variables does not cause serious statistical errors. 

After the initial basis states are generated in the above 
procedure, we have to modify these states to approach 
the true ground state. The ordinary PIRG iteration pro- 
cess, which is denoted as Iteration process B, is based on 
the local update. Therefore, the optimized ground state 
is sometimes trapped in a local minimum of the energy, 
and could be far from the true ground state. To over- 
come this problem, we introduce the following iteration 
process denoted as A, which enables us to modify the 
basis state globally, in spite of the local update. 

The detail of both processes is summarized as follows. 

• Iteration process A 

— To approach the true ground state, we have to 
modify the basis states {|0 a }}- We pick up one 
basis state \4> a ), and compare it with its modified 
state \<f/ a ) = e - AT ^/2 V({s}) e- Arli ^ 2 \cj)a). If 
the energy for \<j>' a ) with appropriate sets of {s} is 
lower than that for \(j) a ), we update the basis state. 
We perform similar updating procedures succes- 
sively for the L basis states. 

• Iteration process B 

— We perform the ordinary PIRG iteration proce- 
dure. 21 ' 22 We first consider a new basis state by 
acting the following three operators on a certain 
basis state \<p a ) step by step: first e~ ATHk ' 2 , an( j 
then e - Ar,7 "*T™.i for all i, and finally e ~ Ar ^/ 2 . 
In each step, we pick up the lowest-energy state 
among the original state and the newly generated 
states. We perform similar procedures for all the 
basis states. 

In our PIRG calculation, we first perform the process A 
a couple dozen times to modify the basis states globally. 
After that, the basis states are locally modified by the 
process B. The combined scheme, which is referred to 
as 'Iteration A — * B', allows us to approach the ground 
state efficiently in small iteration steps and thus obtain 
the physical quantities precisely. 

Furthermore, to improve numerical accuracy, we also 
use the quantum-number projection scheme in the frame- 
work of the PIRG method (PIRG+QP). 25 In this 
paper, we employ a projection operator to the to- 
tal spin-singlet state, Ps=o, which acts on the con- 
verged state obtained in the above procedure: | - 0s=o) = 
Sa=i c ' a Ps=a\4>a) ■ Here c' a is an amplitude reevaluated 
from the generalized eigenvalue problem for a new set of 
the quantum-number projected basis states {Ps=a\4'a)}- 



The PIRG+QP method is particularly powerful to pre- 
cisely determine the phase diagram of frustrated electron 
systems. 26 




PIRG iteration step 

Fig. 3. (Color online) Energy as a function of the number of 
PIRG iteration at L = 500 in the Hubbard model on checker- 
board lattice of N — 16 with U = lOt. Circles, triangles 
and squares represent the results obtained from different initial 
states. For comparison, the results for several different iteration 
processes are shown. Note that the iteration mode is changed 
from A to B between the 9th and 10th step in the case of 'Iter- 
ation A— »B'. 



Here, we demonstrate how effective our method intro- 
duced here is. In Fig. 3, we show the energy as a function 
of the number of iteration for the small cluster of N = 16 
with U = lOt, which is obtained by the PIRG method 
with L = 500, At = 0.5, and At' = At x 10~ 3 . It is 
found that the energy converges in a dozen steps of it- 
eration in spite of the fact that a large number of basis 
states are treated. This is in contrast to the results of 
the ordinary PIRG method where a few hundred steps 
of iteration are needed for the convergence. Note that 
the converged ground state is not necessarily optimized 
by the initial state obtained by the UHF approximation 
with U = Vint - For example, in the present case U = lOt, 
the optimized ground state is deduced more effectively 
starting from the UHF approximation with Ui n t = 4i. 
For comparison, we also show the results obtained only 
by the Iteration B process. It is found that the energy is 
always slightly higher than the energy obtained by 'Itera- 
tion A — > B'. Although the energy difference seems small, 
this optimization process of the ground state plays a cru- 
cial role in deducing the physical quantities precisely. We 
show the results for the cluster of N = 16 with U = lOt 
obtained by the PIRG+QP method with a few choices of 
Uint - In fact, it is found that when the energy variance is 
small, the data in each Ui n t are well fitted to a straight 
line, as shown in Fig. 4. We thus obtain the ground state 
energy and the double occupancy precisely, and confirm 
that the PIRG+QP results are indeed in good agreement 
with those of the exact diagonalization. 

In the following, we carry out PIRG+QP calculations 
for the half-filled systems with N = 16, 20, 24, 28, and 
32 sites in periodic boundary conditions (see Fig. 1). We 
note that in the non-interacting case, the cluster with 
N = 20 has a closed shell structure while the others an 
open shell structure. Therefore, in the latter, the highest 
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Fig. 4. (Color online) (a) Energy and (b) double occupancy per 
site as a function of the energy variance with various initial states 
in the Hubbard model on checkerboard lattice (JV = 16). For 
comparison, exact results are shown as solid symbols. 



energy level in the lower band and the energy level of 
the flat band are degenerate. This means that the Fermi 
point appears only at k = kp[— (jr, 7r)] in the normal 
metallic state. Although this " quasi- Fermi" point arti- 
ficially appears under the periodic boundary condition, 
it might be a good probe to observe the MI transition. 
In the PIRG method, we use At x U/t = 0.5 for both 
iteration processes and At' = At x 1CP 3 for generating 
initial states, and repeat the iteration until the energy 
converges under the truncated Hilbert space. We keep 
the Slater basis states up to L = 500. Note that numer- 
ical errors for the finite cluster (thermodynamic limit), 
which will be shown in the next section, mostly arise from 
the energy variance (system size) extrapolation scheme. 

3. Results 

In this section, we present the numerical results for 
the half-filled Hubbard model on the checkerboard lat- 
tice. First, we compute the ground state energy, the dou- 
ble occupancy and the charge excitation gap. The results 
are shown in Fig. 5. We also perform a finite-size scal- 
ing to obtain the results in the thermodynamic limit, 
as shown in the inset, where a scaling form of N~ 3 ^ 2 
is used for the physical quantities. In the nonintcracting 
case, the dispersive band (a = 2) is fully occupied, where 
Eg/N = -It and {dE g /dU)/N = 0.25. The introduc- 
tion of the Coulomb interaction monotonically decreases 
the double occupancy, as expected for the paramagnetic 
metallic state which persists up to U/t ~ 6. 17,18 Further 
increase in the interaction yields a cusp singularity in 
the ground state energy and a jump singularity in the 
double occupancy, as shown in Figs. 5 (a) and (b). This 
suggests the existence of the first-order phase transition. 
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Fig. 5. (Color online) (a) The ground state energy E g /N, (b) the 
double occupancy (dE g /dU)/N and (c) the charge excitation 
gap A c as a function of U/t on N = 16 (circle) and 32 (triangle) 
lattice and in the thermodynamic limit (square) . Thick gray lines 
in (a) and (b) are the results obtained by the weak and strong 
coupling techniques (see text). Just N = 16 sample is obtained 
by exact diagonalization. N = oo results are obtained from finite 
size scaling, some of which are shown in the inset for U/t = 4 
(circle), 6 (triangle), 8 (square) and 10 (diamond). We note that 
N = 20 system has closed shell structure at half filling so A c = 
0.382 even U = 0. 



The transition point is estimated as U c /t = 6.75 ± 0.25. 
To discuss the nature of the transition, we also calcu- 
late the charge excitation gap A c , which is defined by a 
difference between two chemical potentials, 



A,= 



2 



(7) 



where M + = [E g {M^ + 1, M i + 1) - E g (M h M 4 )] /2, 
M _ = [EgiM^Md-EgiM,- l,Mi-l)]/2 and 
Eg (Mf , Mi ) is the ground state energy of the system, 
where M a is the number of electrons with spin a. We 
find in Fig. 5 (c) that the charge gap is zero for U < U c 
while it is finite for U > U c , i.e. the system is driven to an 
insulating state. It is, however, difficult to quantitatively 
estimate the charge gap in the thermodynamic limit from 
the finite-size scaling analysis within the present numer- 
ical accuracy, as seen in the inset of Fig. 5 (c). 

We now discuss the instability of the insulating state 
against conventional spin or charge ordered states. To 
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this end, we examine the spin and charge correlations 
for the system with the largest cluster of N = 32, where 
the first-order transition occurs near U c (N — > oo) in the 
thermodynamic limit. We first calculate the momentum 
distribution function defined by 



i a (k) = 



2N Yj^kaa^kaa) for fc f (tT.Tt), 

1 " (8) 
Tjvr E^Ux^) for k = ( 7r ' ?r )- 



cr,/3 



The obtained results are shown in Fig. 6. Since the non- 
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Fig. 6. (Color online) The momentum distribution function 
n a (k) in the system with N = 32 for several choices of U/t 
along the symmetry lines in the B.z. 



interacting system has the quasi- Fermi point at (7r, 7r), a 
discontinuity should appear in the momentum distribu- 
tion for the dispersive band n2{k). We indeed observe it 
for U < U c , while this singularity disappears for U > U c , 
in accordance with the fact that the MI transition occurs 
at U = U c , as discussed above. 

We next consider the equal-time correlation functions 
in spin and charge sectors, which are given by 



s mm ' (<?) 



2 

3N 



N/2 



; ' & jm' ) 



N mm >(q) 



x e 



iO -{Rim-R 

„ N/2 

2 / \ 

q-(R im -R , ) 



(9) 



(10) 



where fiim = ni m \ + hi m i and Ri m represents the po- 
sition of the i-th unit cell in the m-th sublattice. Diag- 
onalizing the 2x2 matrix, we obtain S a (q) and N a (q) 
(a = max, min), as shown in Fig. 7. It is found that 
spin and charge correlations are little changed in the 
case U < U c . On the other hand, further increase in the 
interaction U leads to totally different behavior in the 
spin correlation function. It is found that when U > U c , 



Smax(q) is enhanced at q = (0,0) although it never di- 
verges even in the thermodynamic limit. This may sug- 
gest that short-range spin correlations are enhanced in 
the insulating phase. To clarify this point, we also calcu- 
late the site-dependent spin correlation function defined 

by 



C s (n) = 



1 1 



N N n 

i = l T„ = l 



Si+r n ), (11) 



where r„ labels the n-th neighbor site connected by 
transfer integral t and N n is the number of them (A^i = 4, 
N 2 = 2). We find that both C s (l) and C s (2) are always 
negative as shown in the inset of Fig. 7 (a) . Note that the 
spin correlations C s {2) are negative even in the strong 
coupling limit, implying that antifcrromagnctic correla- 
tions are suppressed due to geometrical frustration. As 




(0,0) (0,71) 



(%,%) (71,0) 



(0,0) (7C.7C) 



Fig. 7. (Color online) (a) [(b)] Equal-time correlation function 
S a (q) [N a (q)] in the system with N = 32 for different choices of 
U/t along the symmetry lines in the B.z. The inset of (a) shows 
spin correlation C 3 (n). iV max (0, 0) as a function of U/t is shown 
in the inset of (b). We note that 4S a (q) = N a (q) at U/t = 0. 



for the charge sector, the increase in U/t gradually sup- 
presses charge correlations, as shown in Fig. 7 (b). In 
addition, we find the monotonic decrease in 7V max (0, 0) 
in the inset. Therefore, there is no tendency to stabilize 
the charge ordered state with q = (0,0). 6 

To discuss the spin configuration of the insulating 
phase in detail, we calculate the dimer correlation func- 
tion D a p and the plaquette correlation function Pa(b) 
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Fig. 8. (Color online) (a)-(c) represent the dimer ordering pat- 
terns. Thick (double) lines represent positive (negative) signs for 
the dimer correlation function eq. (12). (d) and (e) represent 
the plaquette ordering patterns. Filled (open) circles represent 
positive (negative) signs for the plaquette correlation function 
eq. (14). (f) represents the n-th neighbor plaquettes, where the 
reference plaquette is labeled by zero. 

defined by 

D a p = {O a 6p), (12) 
1 N 

O a = - for a = 0(13) 

i=l 

pattern a 

and 

Pa { b) = (Q 2 a {B )), (14) 
2 N/2 

Qa [b) = ^ E c- 1 )^. ( 15 ) 

i=l 

pattern A(B) 

where di{= Si ■ Si+$) is the i-th dimer operator and pi[= 
(S ai + S 7i ) ■ (S / 3 i + Ssi)} is the iih plaquette operator. 
The patterns for dimers and plaquettes, and their signs 
(— I) 1 are schematically shown in Figs. 8 (a)-(e). 

Figure 9 (a) shows the U/t dependence of the dimer 
correlations for several possible configurations. We find 
that the dimer correlations are negligibly small for all 
the patterns when U < U c . On the other hand, beyond 
U = U c , the dimer correlations for some patterns are 
suddenly enhanced, while remains small. Namely, 
a crossed-dimcr valence-bond ordered state [see Fig. 8 
(c)], which is realized in the weakly coupled Hcisenbcrg 
chains on the checkerboard lattice, 27 is not stabilized in 
the strong coupling limit. The existence of three dimer 
correlations concludes that the ordinary dimer ordered 
state is not stabilized but the plaquette ordered state 
with pattern A emerges instead. It is indeed seen that 
the positive D^ v appears in the plaquette ordered state 
with A configuration. 

To confirm the above results, we also calculate the 
plaqucttc-singlet correlations, as shown in Fig. 9 (b). It 
is seen that the plaquette correlations with pattern A 



0.02 




Fig. 9. (Color online) (a) The dimer correlation functions D a p 
and (b) the plaquette correlation functions Pa(B) on the N = 32 
lattice as a function of U/t. The inset in (b) is the n-th neighbor 
plaquette correlation C p (n). 



arc enhanced in the region U > U c . Therefore, the quan- 
tum phase transition breaks the translational symmetry, 
leading to the formation of a plaqucttc-singlet ordered 
state. In fact, the plaquette order parameter alternates 
spatially in the insulating phase (U > U c ), as shown in 
the inset of Fig. 9. Here, the n-th neighbor plaquette 
correlations are defined by 

2 1 N/2 N " 

C p( n ) = ^jEE {PiPi+rJ, ( 16 ) 
i— 1 r n — l 

where r„ runs the n-th neighbor plaquettes and N n is the 
number of them (Fig. 8 (f)). These results are consistent 
with those in the Heisenberg limit, 19 where the plaquette 
valence-bond crystal phase is stabilized. We thus end up 
with the conclusion that in the Hubbard model on the 
checkerboard lattice, the normal metallic state is realized 
for small U, while the plaquette-ordered insulating state 
is for large U. The first-order transition occurs between 
these phases at U c /t = 6.75 ± 0.25. 

Here, we wish to make a brief comment on our extrapo- 
lation scheme to the thermodynamic limit. In this paper, 
we have carried out the extrapolation in a scaling form of 
N~ 3 / 2 for the physical quantities, which yields accurate 
results. In fact, we find in Fig. 5 that the obtained results 
are in good agreement with those obtained by the weak 
and strong coupling techniques such as the Green's func- 
tion approach with the self-energy up to second-order in 
U and the fourth-order plaquette expansion around the 
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configuration shown in Fig. 8 (d). We thus confirm that 
our PIRG method with a new iteration scheme works 
well to obtain the reliable results for the Hubbard model 
on the checkerboard lattice. 

4. Summary 

We have studied the Hubbard model on the checker- 
board lattice at half filling by means of the PIRG 
method. When the method has been naively applied 
to our frustrated model, we have encountered a serious 
problem that the approximate ground state is sometimes 
trapped in local minima, leading to unreliable results. To 
overcome this, we have proposed a new iteration scheme 
in the PIRG method, and have demonstrated that the 
satisfactory convergence is achieved in much smaller it- 
eration steps. We believe that our iteration scheme can 
be generally applied to other frustrated electron systems 
for accelerating the convergence and improving the nu- 
merical accuracy. 

It has been found that the increase in the Hubbard 
interaction yields the first-order metal-insulator transi- 
tion at U c /t = 6.75 ± 0.25, where the jump singularity 
appears in the double occupancy. Furthermore, we have 
calculated the dimer and plaquette correlation functions 
to clarify that the system is driven to the plaquette or- 
dered state. 

We have focused on the isotropic checkerboard lat- 
tice model in this paper. If the magnitude of the di- 
agonal hopping is varied, magnetic correlations are en- 
hanced, which should yield a rich phase diagram. This is- 
sue on the anisotropic checkerboard lattice is important 
to study the role of geometrical frustration systemati- 
cally, which is now under consideration. It also remains 
interesting to investigate finite-temperature Mott tran- 
sitions, which should provide further invaluable infor- 
mation about frustrated electrons on the checkerboard 
lattice. 

Acknowledgment 

The authors thank Y. Imai, and T. Ohashi for valu- 
able discussions. Parts of computations were done at 
the Supercomputer Center at the Institute for Solid 
State Physics, University of Tokyo. The work is partly 



supported by a Grant-in- Aid from the Ministry of Ed- 
ucation, Culture, Sports, Science, and Technology of 
Japan [20740194 (A.K.), 20029013 (N.K.) and 19014013 
(N.K.)]. T.Y. is supported by the Japan Society for the 
Promotion of Science. 



1) S. Kondo, D. C. Johnston, C. A. Swenson, F. Borsa, A. V. 
Mahajan, L. L. Miller, T. Gu, A. I. Goldman, M. B. Maple, 
D. A. Gajcwski, E. J. Freeman, N. R. Dilley, R. P. Dickey, J. 
Merrin, K. Kojima, G. M. Luke, Y. J. Uemura, O. Chmaissem, 
and J. D. Jorgensen: Phys. Rev. Lett. 78 (1997) 3729 . 
T. Takeda, M. Nagata, H. Kobayashi, R. Kanno, Y. Kawamoto, 
M. Takano, T. Kamiyama, F. Izumi and A. W. Sleight: J. Solid 
State Chem. 140 (1998) 182. 

H. Sakai, M. Kato, K. Yoshimura and K. Kosuge: J. Phys. Soc. 
Jpn. 71 (2002) 422. 

4) M. Isoda and S. Mori: J. Phys. Soc. Jpn. 69 (2000) 1509. 

5) S. Fujimoto: Phys. Rev. B 64 (2001) 085102. 
S. Fujimoto: Phys. Rev. Lett. 89 (2002) 226402. 



7 
8 
9 
10 

11 

12 

13 
14 
15 
16 

17 

18 

19 

20 
21 
22 
23 
21 
25 
26 
27 



D. Poilblanc: Phys. Rev. Lett. 93 (2004) 197204. 

C. Lacroix: Can. J. Phys. 79 (2001) 1469. 

N. Shannon: Eur. Phys. J. B 27 (2001) 527. 

P. Fulde, A. N. Yarcsko, A. A. Zvyagin and Y. Grin: Europhys. 

Lett. 54 (2001) 779. 

S. Burdin, D. R. Grempel and A. Georges: Phys. Rev. B 66 
(2002) 045111. 

J. Hopkinson and P. Coleman: Phys. Rev. Lett. 89 (2002) 
267201. 

S. Fujimoto: Phys. Rev. B 65 (2002) 155108. 

H. Tsunetsugu: J. Phys. Soc. Jpn 71 (2002) 1844. 

Y. Yamashita and K. Ueda: Phys. Rev. B 67 (2003) 195107. 

M. S. Laad, L. Craco and E. Miiller-Hartmann: Phys. Rev. B 

67 (2003) 033105. 

T. Yoshioka, A. Koga and N. Kawakami: Physica B 378 (2006) 
294. 

A. Koga, T. Yoshioka, N. Kawakami, and H. Yokoyama: Phys- 
ica C 460 (2007) 1070. 

J.-B. Fouet, M. Mambrini, P. Sindzingre and C. Lhuillier: Phys. 

Rev. B 67 (2003) 054411. 

A. Mielke: J. Phys. A 24 (1991) L73. 

M. Imada and T. Kashima: J. Phys. Soc. Jpn. 69 (2000) 2723. 
T. Kashima and M. Imada: J. Phys. Soc. Jpn. 70 (2001) 2287. 
N. Furukawa and M. Imada: J. Phys. Soc. Jpn. 60 (1991) 3669. 
S. Sorella: Phys. Rev. B 64 (2001) 024512. 
T. Mizusaki and M. Imada: Phys. Rev. B 69 (2004) 125110. 
T. Mizusaki and M. Imada: Phys. Rev. B 74 (2006) 014421. 
O. A. Starykh, A. Furusaki and L. Balents: Phys. Rev. B 72 
(2005) 094416 



